﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace PDE
{
    public class PdeTest
    {
        public double[] h;
        public double[] u;
        public int maxStep;
        public double eps;
        public double[] s;
        public double[] c1;
        public double[] c2;

        public PdeTest()
        {
            maxStep = 0;
            eps = 0;
            h = null;
            u = null;
            s = null;
            c1 = null;
            c2 = null;
        }
        
        public PdeTest(double h0, double eps,double u0, int maxStep)
        {
            this.h = new double[maxStep+1];
            u = new double[maxStep+1];
            this.h[0] = h0;
            u[0] = u0;
            this.eps = eps;
            this.maxStep = maxStep;
            s = new double[maxStep];
            c1 = new double[maxStep + 2];
            c2 = new double[maxStep + 2];
        }
        public double reu(double t)
        {
            return u[0]*Math.Exp(3*t);
        }
        public double fTest(double u)
        {
            return 3*u;
        }

        public double RKtest(double hn, double uSt)
        {
            double uNxt=0;
            

            double k1 = fTest(uSt);
            double k2 = fTest(uSt + hn*k1 / 2);
            double k3 = fTest(uSt + hn * k2 / 2);
            double k4 = fTest(uSt + hn * k3);
            uNxt = uSt + hn * (k1 + 2 * k2 + 2 * k3 + k4) / 6;

            return uNxt;
        }
        public void Test()
        {
            int count = 1;
            double[] uTmp =new double[2];
            c1[0] = 0;
            c2[0] = 0;

            while (count <= maxStep)
            {
                u[count] = RKtest(h[count - 1], u[count - 1]);
                uTmp[0] = RKtest(h[count - 1] / 2, u[count - 1]);
                uTmp[1] = RKtest(h[count - 1] / 2, uTmp[0]);
                s[count - 1] = (uTmp[1] - u[count]) / (Math.Pow(2.0, 4) - 1);

                c1[count] = c1[count - 1];
                c2[count] = c2[count - 1];

                if (Math.Abs(s[count - 1]) < eps / Math.Pow(2.0, 5))
                {
                    h[count] = 2 * h[count - 1];
                    c2[count] = c2[count]+1;
                    count++;                         
                }
                else
                {
                    if (Math.Abs(s[count - 1]) > eps)
                    {
                        h[count - 1] = h[count - 1] / 2;
                        c1[count-1] = c1[count-1] + 1;
                    }
                    else
                    {
                        h[count] = h[count - 1];
                        count++;
                    }
                }
            }

        }


    }

    public class PdeTask1
    {
        public double[] h;
        public double[] u;
        public double[] x;
        public int maxStep;
        public double eps;
        public double[] s;
        public double[] c1;
        public double[] c2;

        public PdeTask1()
        {
            maxStep = 0;
            eps = 0;
            h = null;
            u = null;
            s = null;
            c1 = null;
            c2 = null;
            x = null;
        }

        public PdeTask1(double h0, double eps,double u0, int maxStep, double x0)
        {
            this.h = new double[maxStep+1];
            u = new double[maxStep+1];
            this.h[0] = h0;
            u[0] = u0;
            this.eps = eps;
            this.maxStep = maxStep;
            s = new double[maxStep];
            c1 = new double[maxStep + 2];
            c2 = new double[maxStep + 2];
            x = new double[maxStep + 1];
            x[0] = x0;
        }

        public double fTest(double x, double u)
        {
            return u*u/(2*x + x*x) + u - u*u*Math.Sin(10*x);
        }

        public double RK(double hn, double xSt, double uSt)
        {
            double uNxt = 0;


            double k1 = fTest(xSt,uSt);
            double k2 = fTest(xSt+ hn/2, uSt + hn * k1 / 2);
            double k3 = fTest(xSt+ hn/2, uSt + hn * k2 / 2);
            double k4 = fTest(xSt+ hn, uSt + hn * k3);
            uNxt = uSt + hn * (k1 + 2 * k2 + 2 * k3 + k4) / 6;

            return uNxt;
        }

        public void Task1()
        {
            int count = 1;
            double[] uTmp = new double[2];
            c1[0] = 0;
            c2[0] = 0;

            while (count <= maxStep)
            {
                u[count] = RK(h[count - 1],x[count-1], u[count - 1]);
                uTmp[0] = RK(h[count - 1] / 2,x[count-1], u[count - 1]);
                uTmp[1] = RK(h[count - 1] / 2,x[count-1] + h[count - 1] / 2, uTmp[0]);
                s[count - 1] = (uTmp[1] - u[count]) / (Math.Pow(2.0, 4) - 1);

                c1[count] = c1[count - 1];
                c2[count] = c2[count - 1];

                if (Math.Abs(s[count - 1]) < eps / Math.Pow(2.0, 5))
                {
                    h[count] = 2 * h[count - 1];
                    x[count] = x[count -1] + h[count - 1];
                    c2[count] = c2[count] + 1;
                    count++;
                }
                else
                {
                    if (Math.Abs(s[count - 1]) > eps)
                    {
                        h[count - 1] = h[count - 1] / 2;
                        c1[count - 1] = c1[count - 1] + 1;
                    }
                    else
                    {
                        h[count] = h[count - 1];
                        x[count] = x[count -1] + h[count - 1];
                        count++;
                    }
                }
            }
        }


    }


    public class PdeTask2
    {
        public double[] h;
        public double[,] u;
        public double[] x;
        public int maxStep;
        public double eps;
        public double[] s;
        public double[] c1;
        public double[] c2;

        public PdeTask2()
        {
            maxStep = 0;
            eps = 0;
            h = null;
            u = null;
            s = null;
            c1 = null;
            c2 = null;
            x = null;
        }

        public PdeTask2(double h0, double eps,double u0,double v0, int maxStep, double x0)
        {
            this.h = new double[maxStep+1];
            u = new double[maxStep+1,2];
            this.h[0] = h0;
            u[0,0] = u0;
            u[0, 1] = v0;
            this.eps = eps;
            this.maxStep = maxStep;
            s = new double[maxStep];
            c1 = new double[maxStep + 2];
            c2 = new double[maxStep + 2];
            x = new double[maxStep + 1];
            x[0] = x0;
        }

        public double[] fTest(double x, double u, double v, double a) //u`, v`
        {
            double[] tmp = new double[2];
            tmp[0] = v;
            tmp[1] = a * Math.Sqrt(v * v + 1);
            return tmp;
        }

        public double[] RK(double hn, double xSt, double uSt, double vSt, double a)
        {
            double[] uNxt = new double[2];
            

            double[] k1 = new double[2];
            k1[0] = fTest(xSt, uSt, vSt, a)[0];
            k1[1] = fTest(xSt, uSt, vSt, a)[1];
            double[] k2 = new double[2];
            k2[0] = fTest(xSt + hn / 2, uSt + hn * k1[0] / 2, vSt + hn* k1[1]/2,a)[0];
            k2[1] = fTest(xSt + hn / 2, uSt + hn * k1[0] / 2, vSt + hn * k1[1]/2, a)[1];
            double[] k3 = new double[2];
            k3[0] = fTest(xSt + hn / 2, uSt + hn * k2[0] / 2, vSt + hn* k2[1]/2 , a)[0];
            k3[1] = fTest(xSt + hn / 2, uSt + hn * k2[0] / 2, vSt + hn * k2[1] / 2, a)[1];
            double[] k4 = new double[2];
            k4[0] = fTest(xSt + hn, uSt + hn * k3[0], vSt + hn*k3[1], a)[0];
            k4[1] = fTest(xSt + hn, uSt + hn * k3[0], vSt + hn * k3[1], a)[1];
            uNxt[0] = uSt + hn * (k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6;
            uNxt[1] = vSt + hn * (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6;

            return uNxt;
        }

        public void Task2(double a)
        {
            int count = 1;
            double[,] uTmp = new double[2,2]; //номер половинного шага, номер переменной (u - первая)
            c1[0] = 0;
            c2[0] = 0;

            while (count <= maxStep)
            {
                u[count,0] = RK(h[count - 1], x[count - 1], u[count - 1,0], u[count-1,1],a)[0];
                u[count, 1] = RK(h[count - 1], x[count - 1], u[count - 1, 0], u[count - 1, 1], a)[1];
                uTmp[0,0] = RK(h[count - 1] / 2, x[count - 1], u[count - 1,0],u[count - 1,1],a)[0];
                uTmp[0, 1] = RK(h[count - 1] / 2, x[count - 1], u[count - 1, 0], u[count - 1, 1], a)[1];
                uTmp[1,0] = RK(h[count - 1] / 2, x[count - 1] + h[count - 1] / 2, uTmp[0,0], uTmp[0,1],a)[0];
                uTmp[1, 1] = RK(h[count - 1] / 2, x[count - 1] + h[count - 1] / 2, uTmp[0, 0], uTmp[0, 1], a)[1];
                s[count - 1] = (uTmp[1,0] - u[count,0]) / (Math.Pow(2.0, 4) - 1);
                if(Math.Abs(uTmp[1,0] - u[count,0])<Math.Abs(uTmp[1,1] - u[count,1]))
                    s[count - 1] = (uTmp[1, 1] - u[count, 1]) / (Math.Pow(2.0, 4) - 1);

                c1[count] = c1[count - 1];
                c2[count] = c2[count - 1];

                if (Math.Abs(s[count - 1]) < eps / Math.Pow(2.0, 5))
                {
                    h[count] = 2 * h[count - 1];
                    x[count] = x[count - 1] + h[count - 1];
                    c2[count] = c2[count] + 1;
                    count++;
                }
                else
                {
                    if (Math.Abs(s[count - 1]) > eps)
                    {
                        h[count - 1] = h[count - 1] / 2;
                        c1[count - 1] = c1[count - 1] + 1;
                    }
                    else
                    {
                        h[count] = h[count - 1];
                        x[count] = x[count - 1] + h[count - 1];
                        count++;
                    }
                }
            }
        }
    }
}
